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摘 要 : 针对 基于 GPU 求解 大 规模 稀疏 线性 方程 组 进行 了 研究 ， 提 出 一 种 稀疏 和 珑 阵 的 分 块 存储 格式 HMEC (hybrid 
multiple ELL and CSR)。 通 过 重 排序 优化 系数 矩阵 的 存储 结构 ， 将 系数 矩阵 以 一 定 的 比例 分 块 存储 ， 采 用 ELL 与 CSR 
存储 格式 相 结合 的 方式 以 适应 不 同 的 分 块 特征 , 分 别 使 用 适用 于 不 对 称 和 矩阵 的 不 完全 LU 分 解 预 处 理 BICGStab 法 和 对 
称 正 定 和 焉 阵 的 不 完全 Cholesky 分 解 预 处 理 共 斩 梯度 法 求解 大 规模 稀 芯 线性 系统 。 实 验 表 明 ， 应 用 HMEC 格式 存储 稀 
JL AE EH VAIA A GPU kernel 的 方式 实现 前 述 两 种 方法 , 与 其 他 存储 格式 的 实现 方式 作 比 较 , 最 优 可 分 别 获得 31.899% 和 
17.509650 Ja3& 3k XX 
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Optimization of solving sparse linear system based on blocked storage format 


Cheng Kai, Tian Jin, Wu Fei, Wang Ru, Li Hongqin 
(College of Electronic & Electrical Engineering, Shanghai University of Engineering Science, Shanghai 201620, China) 


Abstract: This paper proposed a storage format HMEC (Hybrid Multiple ELL and CSR) of sparse matrix to solve large sparse 


linear equations on GPU. Firstly, we optimized the storage structure of the coefficient matrix by reordering. Secondly, we 
stored the coefficient matrix in a certain scale block. Then we adopt an approach by combining ELL and CSR storage format 


to adapt to different characteristics of blocks. At last, we took Bi-Conjugate Gradient Stabilized (BICGStab) and Conjugate 


p i Gradient (CG) iterative methods to solve large sparse linear systems, they are respectively preconditioned by incomplete-LU 


] s - and incomplete-Cholesky factorization for asymmetric and symmetric positive definite linear matrices. Experiments show that 


comparing the way by storing sparse matrices in HMEC format with other ways by storing in the common storage format, the 
acceleration of the best available we can get are 31.89% and 17.50%. 
Key words: GPU acceleration; conjugate gradient; bi-conjugate gradient stabilized; reorder; HMEC (hybrid multiple ELL and 


CSR) storage format; sparse matrix-vector multiplication 


和 Cusparsel12l 库 的 迭代 法 求解 大 型 方程 组 ， ceu 3 E HH — 
种 基于 HYBI141(Hybrid ) 格 式 的 混合 存储 格式 ， 通 过 持续 划分 
求解 大 规模 稀 疏 线性 系统 广泛 应 用 于 计算 科学 与 工程 的 各 ELL 以 提升 性 能 ; 阳 王 东 等 人 针对 SPMV 的 优化 提出 了 一 种 重 
个 领域 0421。 目 前 ， 常 用 的 求解 稀 疏 线性 系统 的 方法 有 直接 法 排序 的 分 块 存储 格式 BCEI151， 良 好 地 适应 了 有 不 规则 非 零 元 
与 迭代 法 3.41， 其 中 ， 和 迭代 法 尤为 适合 大 规模 稀 朴 线性 系统 的 结构 的 矩阵。 本 文 在 阳 王 东 的 研究 基础 上 将 系数 矩阵 4 重新 排 
求解 。 常 用 的 迭代 法 包括 固定 的 迭代 法 ， 如 基于 分 割 的 Jacobi FF, 在 殷 建 的 研究 上 采用 多 次 划分 ELL 格式 以 减弱 并 行 计算 时 
和 Gauss-Seidel (GS) 法 ;不 固定 的 迭代 法 ,如 Krylov 子 空间 方法 ， ”的 负载 不 平衡 ， 并 采用 CSR 格式 存储 剩余 矩阵 以 提升 存储 效 
包括 可 适用 于 不 对 称 系统 的 稳定 双 共 轿 梯 度 (bi-conjugate — 率 。 最 终 实验 表明 ， 采 用 本 文 所 提出 的 分 块 存储 格式 对 于 大 规 
gradient stabilized, BICGSTAB) (516117l 法 和 对 称 正定 线性 系统 模 稀疏 线性 系统 的 求解 具有 良好 的 加 速效 果 。 
WHAE (conjugate gradienb) 法 [81。 在 实践 中 ， 当 稀疏 矩阵 为 

病态 矩阵 时 , 需 用 预 处 理 技术 加 以 改进 以 达到 良好 的 计算 效果 。 1 稀 玻 矩阵 的 存储 格式 

Mayer 9! 研究 了 基于 多 级 不 完全 LU 分 解 预 处 理 法 求解 线 因 稀 疏 和 矩阵 中 非 零 元 素 分 布 的 不 规则 性 


limi 


IN 


导致 采用 单一 的 
性 三 角 系 统 ; NaumovI10] 在 Mayer 的 基础 上 研究 运用 Cubalslll! 存储 格式 进行 运算 时 ， 取 得 的 效果 并 不 理想 ， 鉴 于 此 ， 多 种 优 
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| 
录用 稿 程 凯 ， 等 : 基于 分 块 存储 格式 的 稀 路 线 性 系 乡 天 化 


化 矩阵 存储 格式 的 方法 相继 提出 。 一 种 是 对 行 合 并 分 块 以 缓解 ilililalala 111 0| 1 
各 行 间 非 零 元 数 差异 过 大 , 如 分 块 的 CSR 方法 (BCSR) 1617, PUDE E Ut = E add 
(BELLPACK) !18! f& 3X 5$, 4; Hr Ry ELLPACK 格式 : Saco ius I ! ! - : 
(SELL-PACK)I191 等 ;, 另 一 种 是 采取 混合 格式 进行 存储 ， 以 适应 1 1 0 
不 规则 的 稀 朴 特征 ,如 混合 ELLIELLPACKJ 和 COO(Coordinate) 1 1 0 


格式 的 HY BÉERUSEXSE. ] 存 储 格 式 等 。 FEAE RE 值 "— g 
AX E Ua FUR RB EET Be tr fé X HMEC， 由 多 个 
ELL 块 和 一 个 CSR Hii f Ifi JE AKARE AE RE TJ 


vas | 1|[2|]3|4|5/6|7|8 9]|10 


COO 部 分 Row 0|10101011111111 2 2 


多 个 子 和 矩阵 分 别 计算 ， 并 采用 重新 排序 的 方法 改善 稀疏 矩阵 的 
非 零 元 结构 。 其 中 ，CSR 块 存储 分 割 出 的 离散 点 ， 良 好 地 适应 cm |2|[3]4|5]2/3|4|5|2.3 
Y EJBOBEERIASUUTE, BUKBUSESTE f Mio ERIT ERNE . 图 3 HYB 存储 格式 


| 2.3 HMEC 存储 格式 
人 本 文 提出 的 HMEC 存储 格式 由 多 个 ELL 块 和 一 个 CSR X 
如 图 1 所 示 ，ELL 格式 用 两 个 数组 values 和 col 来 存储 ”混合 而 成 ,即将 经 过 重 排序 的 矩阵 4 在 未 达到 停止 划分 的 条 件 
一 个 4x2 RUBER, 其 中 , 数组 values Fi REER RETE ” 前 , 根据 闵 值 持续 划分 为 ELL 阵 和 剩余 矩阵 ， 最 终 将 矩阵 A 
v— o 零 元 的 值 ， 数 组 col 存储 相应 非 零 元 素 在 原 系数 矩阵 中 的 列 索 ”保存 为 多 个 以 ELL 格式 存储 的 矩阵 块 和 一 个 以 CSR 块 格式 存 
” B. 储 的 矩阵 块 ， 每 个 格子 代表 矩阵 的 一 个 元 素 ， 空 格 代表 0。 为 


1 p 1 | 2 0 | 2 便于 简 述 ， 以 划分 2 次 ELL 为 例 ， 如 图 4 所 示 。 
3 3 | 0 1 | 0 eu EuD CSR 
4 4|0 2 |0 | | | 
5 6 NA AES "bi | du 
: EEE Values Col 1 | m 4|1|1| 
= 图 1 ELL 存储 格式 3|11111 | 一 p S 
© 21 CSR 存储 格式 EE mctu 
N 如 图 2 所 示 , CSR 格式 采用 三 个 数组 value, rowptr, colind —| na 
2^ Riik ABSWRRE, IH values (AREE RAE - - 
零 元 的 值 ， 数 组 colind efl T 3:975 BURG Pe rb D A m 
引 ， 数 组 rowptr 存储 了 每 一 行 元 素 在 压缩 存储 格式 中 的 起 始 位 1 1.1/1 1 第 二 次 划分 = 
m. 1|1 o 四 加 
1 2 Values 1:213]4[|5/|6 -—- ea 
3 
" Rowptr 023.416 图 4 HMEC 存储 格式 
5 6 对 于 系数 矩阵 A， 为 优化 矩阵 A 内 的 非 零 元 分 布 ， 采 用 选 
Conad 0[2:11[2:1/|3 择 排 序 法 按 每 行 非 零 元 的 个 数 由 高 至 低 进行 降序 排序 ， 并 额外 
全民 用 一 个 矩阵 recoder 记录 下 排序 时 的 行 索引 ， 可 得 矩阵 41。 将 
图 2 CSR 存储 格式 FERE A1 依据 阔 值 天 划分 为 一 个 ELL 格式 存储 的 矩阵 ELLI 和 


2.0 HYB 存储 格式 剩余 矩阵 42， 将 矩阵 42 fills BUR K RIISI — A ELL 格式 存储 

如 图 3 所 示 , 为 便于 简 述 ， 以 一 个 6 阶 方 阵 4 为 例 , HYB — RAEEE ELL2 和 剩余 矩阵 A3， 划 分 完毕 后 ,矩阵 43 以 CSR 格 
存储 格式 依据 阔 值 天 将 矩阵 划分 为 ELL 格式 与 COO 格式 两 部 。 式 存储 ， 最 终 划 分 为 ELLI. ELL2 和 CSR 三 部 分 并 分 别 以 相 
分 存储 ， 每 个 格子 代表 和 抢 阵 的 一 个 元 素 ， 空 格 代表 0。ELL 格 应 的 格式 存储 ， 阔 值 天 的 划分 规则 如 下 


式 存储 矩阵 4 中 非 零 元 的 值 和 它 所 在 的 列 , COO 格式 存储 划分 求 出 矩阵 S 中 第 一 ubica x《 第 一 次 划分 S 为 41， 
后 剩余 矩阵 的 值 、 行 和 列 ， 并 分 别 用 三 个 数组 values, row 和 第 二 次 划分 S 为 A2... 
colind 表示 。 if(x > 512){ 

for j22: 1: sf 


% j 为 矩阵 S 的 列 ，sf 为 矩阵 S 的 总 列 数 
求 出 前 j 列 的 总 非 零 元 数 s; 
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录用 稿 
求 出 前 j 列 的 元 素 的 总 数量 c ; 
若 s> (2c/3), 则 K=j; 否则 break， 跳 出 循环 ; 
o d BOE. K (矩阵 为 41 Wf. K BC KIT, 矩阵 为 42 时 ， 
K Ht K2...) 
} 


其 中 ， 第 一 列 的 非 零 元 数 x>512 为 可 以 持续 划分 的 条 件 ， 
RRE EE rH HERE TUR B UP T 512 时 会 停止 划分 ， 因 为 此 
时 若 继 续 划分 已 经 不 能 够 再 获得 性 能 的 提升 。 


3 ”方程 组 求解 方法 


3.4. HEA 

共 斩 梯 度 法 的 研究 基础 为 构造 搜索 方向 和 步 长 因子 从 而 
找到 合适 的 下 降 方 向 ， 再 进行 迭代 计算 ， 最 终 找 到 最 优 解 并 给 
出 收敛 性 。 
体 算法 可 描述 为 : 对 于 无 约束 优化 问题 , 给 出 一 个 初始 
值 xo, 算法 迭代 产生 亏 , xz, xXx， 若菜 一 x 是 优化 问题 的 解 ， 
在 第 Kk 次 迭代 时 ， 当 前 迭代 点 为 x， 则 线 搜索 型 的 方法 产生 一 
个 搜索 方向 pr 后 ， 下 一 个 迭代 点 xen. 的 计算 公式 为 : 
Xa =R tP AP a 是 步 长 因子 , 它 满足 某 些 线 搜索 条 件 。 
3.2 不 完全 Cholesky 分 解 

对 于 方程 组 


Ax=b (1) 
其 中 : 4 n 阶 对 称 正定 矩阵 , b 为 列 向 量 ,x 为 待 求解 。 若 对 
HERE A 进行 不 完全 Cholesky 分 解 ， 可 得 式 (2)， 其 中 工 为 下 三 
HERE, D IIAIR, RAEE E IEE A £I LDLT Z 25. 
若 矩 阵 4 中 零 元 素 4; =0， 可 人 为 地 令 和 矩阵 工 中 匀 =0 ， 以 使 
和 矩阵 工 的 稀 跑 性 与 矩阵 A 一致。 
A=LDĽ +E (2) 
EHER LDO ERARO), TRG). 


(LD'^y" A(u.D')')" (LD'?y x - (LD*)p (3) 


令 式 (3) 中 ADY =R, Ht, 
处 理 矩 阵 M 如 式 (4) 所 示 。 


玉 为 上 三 角 和 矩阵 ， 可 得 预 


A~M = R’'R (4) 

式 (3) 可 简化 为 式 (5)。 
(RT AR (Rx) = (R7 b) (5) 
435) (RTAR”) =C, Rx=y, (R7b)=p , 可 得 式 (6)。 
Cy-p (6) 


式 (6) 中 C 即 为 经 过 不 完全 Cholesky 分 解法 预 处 理 的 系数 


ABE. 
3.8 BICGSTAB 算法 

BICGSTAB 算法 的 主要 思想 是 : 对 于 n 阶 线性 方程 组 
Ax=b, 令 初始 近似 解 为 xo, 第 外 次 近似 解 为 xe HARER 
Xon b - Ax, 有 k Br Krylov 子 空间 


K, =span{n, An, A" n), K, 2span(4, A75, (A^) i} 。 该 算 


i fEKrylov-T- 25 [8] KeFP 3e UP Ype 通过 选择 合适 的 参数 w 和 


Bo FERM, Hilt rex * Kr LK, 计算 的 残 


25 总 是 取 最 小 范 数 。 有 : na mrt OP Pen = ha TAP o 
在 求解 实际 问题 中 ，BICGSTAB 算法 通常 与 预 处 理 技术 
相 结 合 以 提升 收敛 速度 ， 增 强 稳定 性 。 流 程 如 下 : 

a) 给 定 系数 矩阵 4， 向 量 b，， 初 始 值 wo， 最 大 迭代 次 数 kiwar 
最 大 容许 误差 Ena 和 预 处 理 和 矩阵 M。 则 初始 残 差 m= -Axo 令 
kel, =h 。 

bD) Kk S kpa HE Eno MEKEEC), 否则 , 停止 并 输出 六 。 
车 Pi=0 或 @.1=0，, 算法 终止 , 输出 失 


c) 4a 2 P nas 
败 信息 ， 和 否则 转 d)。 
d k-LNlp -n, mW i 
Bia = G0 4043 / G9, 3044) Dy = Tia + Piana ~ OM) 
6) 求解 方程 MP = p , VE SEV, = Apo, = p, 1G V), 


Sk = hpa T AV; o 


Ð e=|si|, FEZ Emo Mx -x,-o,P,. B, 停止 输 


E x o 
DREMS =s, t=4A8, =t sitt), 


X=X tA P, +tOS, n—s—-ÓOM, & - |n. 


， 令 k=k+1， 转 


b) 。 
3.4 不 完全 LU 分 解 

运用 不 完全 LU 分 解 获取 式 (1) 的 预 处 理 和 矩阵 M。 若 矩阵 A 
的 顺序 主子 式 都 非 奇异 ， 则 一 定 能 进行 LU 分 解 。 若 取 预 处 理 
矩阵 直接 进行 LU 分 解 ， 得 M=LU， 则 理论 上 最 优 ， 但 在 实际 
应 用 中 并 不 可 行 ， 一 种 不 完全 LU 分 解 是 指 把 系数 矩阵 A 分 解 
为 和 上 U， 它们 的 两 个 因子 分 别 与 4 的 下 、 上 三 角 部 分 有 完全 
相同 的 非 零 元 结构 ， 这 种 分 解 就 是 无 填充 的 不 完全 LU 分 解 ， 
记 为 ILU(0)， 如 式 (7) 所 示 。“0” 代 表 因 子 分 解 是 产生 的 非 零 元 
素 的 个 数 为 0。 


AxM - LÜ 0 
4 CUDA 平台 的 实现 


4.1 基于 HMEC 格式 的 算法 

HMEC 格式 的 数据 结构 由 多 个 ELL 部 分 和 一 个 CSR 部 分 
组 成 ,其 在 CUDA 平 台 上 的 SpMYV 算法 由 若干 个 ELL kernel. 和 
一 个 CSR kernel 组 成 。 如 图 5 所 示 ， 在 计算 过 程 中 ， 先 将 存储 
JI HMEC 格式 的 稀疏 矩阵 4 和 向 量 x 从 host memory 中 传输 到 
global memory, 再 依次 启动 这 些 kernel， 计 算 完成 后 ， 将 结果 从 
global memory 传 回 host memory. 在 整个 算法 执行 过 程 中 , ELL 
所 需 启动 的 线程 数 随 着 剩余 矩阵 行 数 的 减少 而 减少 。 
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录用 稿 程 Zt, F: - 


代数 子 程序 )。 在 使 用 Cublas library 时 ， 应 用 程序 必须 在 GPU 
| 内 存 空间 中 分 配 所 需 的 矩阵 向 量 空间 ， 并 将 其 填充 数据 ， 再 顺 
Global memory 序 调用 所 需 的 Cublas 函数 ， 最 后 将 计算 结果 从 GPU 内 存 空间 
上 传 到 主机 中 。Cusparse library 包含 了 一 系列 处 理 稀 玻 矩阵 的 


Host memory 


基本 的 线性 代数 子 程式 , 是 CUDA 函数 库 的 一 部 分 , 从 C, C++ 
ELL ini "m eoo ELL kernel FE 中 调 Jo 
本 文 的 SpMYV 求解 问题 正 可 以 归结 为 稀疏 矩阵 与 一 个 稠密 
向 量 间 的 操作 ， 运 用 Cusparse library 可 以 较 好 地 实现 。 对 于 大 
QS iE VE ISARBIIRAE,. ASOCUEHH TOS T TOSUBIOBEEITIAS SE 
Laid 全 LU 分 解 预 处 理 BICGSTAB 法 和 对 称 正 定 矩 阵 的 不 完全 
Y Cholesky 4) fé TRAE ERECTUS EYE, ERE, GPU fá vix 
iu dois 代 前 和 每 次 迭代 过 程 中 的 矩阵 与 向 量 ， 向 量 与 向 量 间 的 并 行 计 
图 5 基于 HMEC 格式 的 执行 算法 算 ，CPU 负责 迭代 循环 和 收敛 条 件 的 控制 以 及 标量 相 除 等 操 
对 于 大 型 稀 玻 线性 系统 , 其 系数 矩阵 4 E, H Matlab 41E. 
对 和 矩阵 A 进行 处 理 ， 以 ELL kernel 为 例 ， 可 得 对 应 格式 的 值 、 运用 Cublas 、Cusparse 库 函 数 进行 求解 的 函数 ; 
列 、 用 于 记录 和 矩阵 A 按 每 行 非 零 元 的 个 数 降序 排序 时 产生 的 行 调用 cusparseCreateMatDescr 函数 创建 descr. L, descr_U, descr. A, 
索引 jds row 及 每 次 划分 ELL 时 的 非 零 元 行 数 和 浆 值 K 并 调用 cusparseSetMatIndexBase 与 cusparseSetMatType 函数 进行 
运用 CUDA 平台 编写 的 SPMV 算法 如 下 : 初始 化 。 
D 1) ELL kernel: 调用 cusparseDcsrsv_analysis 函数 生成 info_L,info_U 
e if (row« ELL 格式 的 行 数 ) 计算 初始 残 差 r=b 一 Axo 
Y. for (int i = 0,1,2... BË K) 调用 cusparseDcsrmv 函数 计算 r=Ax0 
对 应 行 、 列 中 元 素 相 乘 ， 得 到 结 调用 cublasDscal 函数 计算 r—-Axo 
根据 行 索引 数组 jds_row 求 得 经 排序 后 的 结果 调用 cublasDaxpy 函数 计算 r= b+r= b — Axo 
2) CSR kernel: 
pe if (row2«CSR 格式 的 非 零 元 行 数 ) 5 ”实验 分 析 
» for (int elem —row. ptr[row2]...row ptr[row2--1]) 本 文 的 计算 平台 为 CPU: Intel Core™ i5-6500 93.20 GHz; 
. 对 应 行 、 列 中 元 素 相 乘 ， 得 到 结果 GPU: NVIDIA GeForce GTX1050; 运行 内 存 (RAM) 的 大 小 
根据 行 索引 数组 jds_row 求 得 经 排序 后 的 结 为 8 GB; 系统 是 : Windows7 64 位 ， 实 验 环 境 为 : Visual Studio 
: = = 42 基于 Cublas library、Cusparse library 的 算法 2012, CUDA8.0。 如 表 1 所 示 ， 展 现 了 测试 所 用 的 稀疏 矩阵 ， 
f Cublas library 就 是 在 NVIDIA CUDA 中 实现 Blas( 基 本 线性 MAH LAB EARR A UF Sparse Matrix Collection20! , 


表 1 IAH F AE E 


名 称 行 数 非 零 元 数 平均 每 行 非 零 元 数 
HB/685_bus 685 3,249 4.74 
HB/1138_bus 1138 4,054 3.56 
HB/bespwr07 1612 5,824 3.61 
HB/bespwr10 5300 21,842 4.12 

Nasa/barth 6691 26,439 3.95 
HB/gemat12 4929 33,044 6.7 
HB/besstk23 3134 45,178 14.42 
HB/besstk24 3562 159,910 44.89 
HB/besstk28 4410 219,024 49.67 
HB/besstk25 15439 252,241 16.34 
HB/bcsstk16 4884 290,378 59.45 
HB/besstk17 10974 428,650 39.06 


HB/psmigr 1 3140 543,160 173 
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名 称 行 数 非 零 元 数 平均 每 行 非 零 元 数 
HB/bcsstk33 8738 591,904 67.73 
HB/bcsstk29 13992 619,488 44.277 

Hollinger/g7jac180 53370 641,290 12.01 
Boe/crystk02 13965 968,583 69.35 
Boe/bcsstk36 23052 1,143,140 49.59 

ATandT/twotone 120750 1,206,265 9.99 
HB/bcsstk30 28924 2,043,492 70.65 

Boe/ct20stif 52329 2,600,295 49.69 

Rothberg/3dtube 45330 3,213,618 70.89 
Chen/pkustk04 55590 4,218,660 75.89 

Chen/pkustk12 94653 7,512,317 79.37 

Rothberg/gearbox 153746 9,080,404 59.06 
Boeing/pwtk 217918 11,524,432 52.88 

5.1 SpMV 计算 时 间 分 的 次 数 有 合适 的 阔 值 p， 因 本 文 所 采用 的 矩阵 的 行 数 在 六 百 


对 于 大 型 稀 玻 方程 组 , 采用 调用 GPU kernel 的 方式 计算 一 ”到 二 十 万 的 量 级 间 ， 经 本 文 实验 可 知 , p 取 2 或 3 为 宜 , 这 里 p 
次 SpMV， 其 中 ，HEC 存储 格式 表示 划分 一 次 ELL， 将 矩阵 A W2, 即将 HMEC 格式 划分 2 次 ELL 格式 的 方式 , 计算 时 间 如 
划分 为 ELL 和 CSR 两 部 分 存储 。 因 实验 条 件 的 限制 ， 随 着 划 WK 2 所 示 。 
分 次 数 的 增多 ， 采 用 MATLAB 进行 处 理 的 时 间 也 越 长 ， 故 划 
R2 单 次 SpMV 计算 时 间 ( 单 位 ms) 


名 称 COO CSR ELL HYB HEC HMECQ 次 ) 
HB/1138_bus 0.019 0.024 0.022 0.026 0.027 0.028 
HB/bespwr07 0.027 0.029 0.025 0.029 0.036 0.029 
HB/bespwr10 0.052 0.031 0.039 0.047 0.045 0.041 

Nasa/barth 0.058 0.039 0.046 0.063 0.058 0.048 
HB/gemat12 0.062 0.037 0.082 0.035 0.033 0.033 
HB/besstk28 0.192 0.163 0.167 0.129 0.121 0.112 
HB/psmigr 1 2.212 2.054 3.58 1.26 1.871 1.567 
HB/bcsstk33 0.592 0.471 0.495 0.425 0.395 0.362 
HB/besstk29 3 0.57 0.698 0.412 0.386 0.335 
Hollinger/g7jac180 1.25 0.78 2.61 0.642 0.593 0.568 
Boe/crystk02 .48 0.83 4.53 0.715 0.653 0.614 
Boe/besstk36 1.65 0.99 1.425 0.874 0.71 0.645 
ATandT/twotone 3.127 2.854 7.48 2.536 2.215 1.984 
HB/besstk30 1.92 1.68 2.089 1.398 1.112 1.034 
Boe/ct20stif 2.86 2.52 38.59 221 1.88 1.69 
Rothberg/3dtube 4.35 3.95 31.14 3.66 3.09 2.97 
Chen/pkustk04 3.89 3.16 68.51 3.27 2.62 2.36 
Chen/pkustk12 10.68 8.59 95.63 8.25 7.14 6.52 
Rothberg/gearbox 2.868 0.952 5.04 3.874 3.498 3.287 
Boeing/pwtk 5.914 5.152 12.92 4.967 4.275 3.876 
表 2 可 知 ， 虽 然 ELL 格式 非常 适合 特征 为 矢量 的 结构 ， SEDE REFTARE, MKE- ERB, ACEH S 
但 是 当 稀 朴 矩阵 每 行 的 非 零 元 素 值 的 数目 变化 特别 大 时 ，ELL — AED UBER ROSE RT, 在 本 实验 中 , R k E 30000 


Hf 


格式 就 会 丧失 这 种 高 效 性 。 当 稀疏 矩阵 4 的 非 零 元 数 较 少 ， 即 ”左右 适宜 。 原因 是 CSR 格式 只 需 存 储 和 矩阵 A 中 非 零 元 的 值 、 列 
小 于 一 定 的 阔 值 x 时 ，CSR 存储 格式 所 需 的 时 间 最 少 ， 当 稀 琉 ”和 行 偏 移 ， 当 非 零 元 数 较 少 时 ， 存 储 的 数据 量 也 较 少 ， 有 较 快 
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的 计算 速度 当 非 零 元 数 较 多 时 ， 随 着 存储 的 数据 量 的 增多 ， 

CSR 格式 的 求解 速度 受到 了 一 定 的 制约 ， 而 HMEC 存储 格式 
因 拥 有 ELL 格式 更 善于 并 行 计算 的 优势 ， 且 经 多 次 划分 ELL 
EAE RETER A 中 非 零 元 的 值 和 列 ， 存 储 的 数据 
量 较 少 ， 拥 有 更 快 的 并 行 计算 速度 。 

在 本 次 单 次 SpMV 运算 实验 中 ,对 比 于 单 次 划分 ELL 格式 ， 
随 着 矩阵 A 中 非 零 元 数 及 大 小 的 增加 , HMEC 格式 相 较 于 HEC 
格式 的 优势 愈加 凸显 ， 对 于 矩阵 HB/psmigr 1， 达到 了 16.25% 
的 加 速效 果 ; 对 比 于 现 有 的 通用 存储 格式 (COO、CSR、ELL.、 
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Chi 
R0 z, PPP N 系统 


合作 期 刊 
程 9X 优化 
应 用 广泛 。 如 在 电磁 问题 中 ， 由 有 限 元 单元 法 离散 化 导出 的 线 


性 代数 方程 组 (其 系数 矩阵 记 为 4) 主要 有 以 下 基本 特征 : 
a) JERE A 是 稀疏 的 。 网 格 节点 一 般 是 很 多 的 ， 但 每 一 

点 上 的 离散 化 方程 只 联系 几 个 相 邻 节点 ， 除 与 这 些 节 Red 

未 知 解 的 系数 不 为 零 外 ， 其 余 系数 都 是 0， 所 以 ， 系 数 矩 阵 的 

元 素 大 量 是 0， 这 就 是 系数 矩阵 的 稀 玻 性 。 


DERE A 是 对 称 且 正定 的 。 
co^ BE A 大 都 是 病态 的 。 若 4 的 条 件数 
cond(A) = (Apax | Arin) > 1 (As 和 um 分 别 为 4 的 最 大 特征 值 


HYB)， 对 于 矩阵 Boe/ct20stif， 本 文 提 出 的 稀 疏 矩阵 存储 格式 
HMEC 有 最 好 的 加 速效 果 ， 相 较 于 CSR 格式 达到 了 32.93%, 
而 对 于 
JERE Rothberg/gearbox, CSR 格式 取得 了 最 好 的 加 速效 果 ， 
这 是 由 于 该 矩阵 的 非 零 元 值 都 为 1， 数 据 量 少 ， 更 利于 CSR 
Kernel 的 并 行 计 算 。 
5.2 不 完全 LU 分 解 预 处 理 BICGSTAB 法 

如 图 6 所 示 ，HMEC 采用 划分 2 次 ELL 格式 的 方式 ， 
CSR(Cusparse) 和 HYB(Cusparse) 分 别 代 表 调 用 cusparseDesrmv() 
和 cusparseDhybmvO 函 数 计算 SPMV 的 方式 ， 运 用 不 完全 LU 
分 解 预 处 理 BICGSTAB 法 计算 时 间 如 图 6 所 示 。 
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图 6 BICGSTAB 法 计算 时 间 

图 6 可 知 , 本 文 提出 的 稀疏 矩阵 存储 格式 HMEC 在 进行 
多 次 的 SPMV ZR, ARRIER. FFEA, FERAE 
KE 4 为 HB/bcsstk33 时 ， 该 存储 格式 有 所 有 格式 中 相对 最 短 的 


GPU 并 行 计 算 时 间 ， 为 16.32 ms， 相 较 于 CSR 存储 格式 加 速 
效果 显著 ， 达 到 了 31.89% 。 


矩阵 HB/besstk33 的 非 零 元 数 为 591904， 远 高 于 和 矩阵 
HB/1138 bus 的 4054， 计算 时 间 却 相近 ， 这 是 由 于 和 矩阵 
HB/bcsstk33 中 的 非 零 元 值 都 为 整数 1, 大 大 加 快 了 SpMYV 运算 
的 求解 速率 。 

和 矩阵 HB/bcesstk16 中 每 行 的 非 零 元 数 都 相同 ， 不 存在 剩余 
ERE, AIMERA ELL 存储 格式 ， 无 法 进行 HYB 等 混合 格式 


与 最 小 特征 值 的 绝对 值 ) 时 ， 称 4 为 病态 矩阵 。 当 A DASE 


阵 时 ， 需 用 多 种 方法 加 速 对 应 方程 组 的 求解 。 
HMEC 采用 划分 2 次 ELL 格式 的 方式 ， 采 用 ICCG 法 的 
计算 时 间 如 图 7 所 示 。 
70 
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E 7 ICCG 法 计算 时 间 
图 7 可 知 ， 本 文 提 出 的 稀 疏 矩阵 存储 格式 HMEC 在 
用 ICCG 法 求解 对 称 正定 线性 系统 中 有 着 良好 的 加 速效 果 。 
WARE A 为 HB/sherman3 时 , HMEC 格式 在 相同 矩阵 中 耗 时 
Hd. 7427.85 ms， 相 较 于 CSR 存储 格式 达到 了 17.50% 的 加 
速 性 能 ， 相 较 于 HEC 格式 达到 了 13.29% 的 加 速效 果 。 


lK bi 
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a) 本 文 提出 的 稀疏 矩阵 混合 存储 格式 HMEC 综合 了 ELL 
和 CSR 存储 格式 各 自 的 优势 ， 通 过 不 断 的 划分 ELL 块 来 降低 
于 计算 的 矩阵 4 的 大 小 , 并 利用 重 排序 改善 矩阵 4 的 非 零 元 
结构 ， 拥 有 良好 的 适应 性 。 

b)HMEC 存储 格式 在 SpMV 计算 中 加 速效 果 显 著 , 尤其 
矩阵 4 较 大 、 非 零 元 数 较 多 时 ， 拥 有 优秀 的 加 速 性 能 ， 达 到 了 
32.9396, 

cy TOS Biz E S ESUOR AE, 2d 
矩阵 时 ,采用 不 完全 LU 分 解 预 处 理 BICGSTAB 法 求解 ， 稳 定 
性 较 高 , 并 达到 了 31.89% 的 加 速效 果 ; H ROB E 4 为 对 称 正 


S 


M 


系数 和 矩阵 4 为 不 对 称 


的 划分 ， 表 中 NA 代表 其 数值 与 ELL 格式 的 时 间 相 同 ， 这 里 

ELL 格式 的 计算 时 间 与 HYB(Kernel) 的 计算 时 间 相 一 致 。 

5.8 不 完全 Cholesky 分 解 预 处 理 共 斩 梯 度 法 
目前 ,不 完全 Cholesky 4 f TRUE IEH S |j RETE (incomplete 


Cholesky factorization preconditioned conjugate gradient, ICCG) 


TEEME, RA ICCG HITRE TEIXERET A E, d 
有 良好 的 加 速效 果 ， 达 到 了 17.50%. 

本 文 从 优化 稀 政 和 矩阵 的 存储 格式 角度 来 加 速 大 型 稀 琉 线 
性 系统 的 求解 ， 但 进行 SPMV 运算 时 ， 采 用 了 单线 程 计 算 和 矩阵 
A 中 一 行 元 素 与 向 量 x 中 一 个 元 素 的 乘积 ，HMEC kernel 函数 
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仍 有 进一步 优化 的 空间 , 如 采用 以 Half-Warp 为 单位 进行 循环 
计算 的 优化 方法 ， 以 减少 线程 的 空转 等 待 ， 或 使 用 纹理 存储 器 
来 存储 稀 玻 矩阵 的 值 ， 以 加 速 数据 读 取 。 此 外 ， 本 文 提 及 了 
ICCG 法 在 电磁 问题 分 析 中 的 应 用 ， 如 何 使 本 文 提出 的 加 速 求 
解 方法 成 功 地 应 用 于 电磁 问题 并 扩大 其 适用 范围 显得 尤为 如 
要 ， 这 将 是 下 一 步 的 研究 方向 。 
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